/*
    Copyright 2008-2011 Patrik Jonsson, sunrise@familjenjonsson.org

    This file is part of Sunrise.

    Sunrise is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 3 of the License, or
    (at your option) any later version.

    Sunrise is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with Sunrise.  If not, see <http://www.gnu.org/licenses/>.

*/

/** \file

    Contains the p04_grain_model, which models the grain model used by
    the Pascucci et al 04 benchmark. */

// $Id$

#ifndef __p04_grain_model__
#define __p04_grain_model__

#include <vector>
#include "grain_collection_grain_model.h"
#include "grain_size.h"
#include "thermal_equilibrium_grain.h"
#include "misc.h"

namespace mcrx {
  template <template<class> class, typename> class p04_grain_model;
  template <template<class> class chromatic_policy, typename T_rng_policy> 
  class p04_grain_model_creator_functor;
};


/** This class implements the grain model used by the Pascucci et al
    04 benchmark. It is a straightforward derivation of
    grain_collection_grain_model which just loads the p04 dust opacity
    file. The class methods are not reentrant and can not be used by
    multiple threads. For the temperature calculation, you must make a
    copy of the object for each thread.  */
template <template<class> class chromatic_policy, 
	  typename T_rng_policy = mcrx::global_random> 
class mcrx::p04_grain_model : 
  public mcrx::grain_collection_grain_model<chromatic_policy, T_rng_policy> {

  void retarded (const std::string& data_directory, const Preferences& p,
		 const T_unit_map& units);

public:
  p04_grain_model (const std::string&, const Preferences& p, 
		   const T_unit_map&);
  p04_grain_model (const p04_grain_model& g) :
    // need to make sure array members are copied. We make the copies
    // thread-local, since the methods are not reentrant anyway.
    grain_collection_grain_model<chromatic_policy, T_rng_policy>(g) {};

  virtual boost::shared_ptr<grain_model<chromatic_policy, T_rng_policy> > 
  clone() const {
    return boost::shared_ptr<grain_model<chromatic_policy, T_rng_policy> >
      (new p04_grain_model(*this));};
};


/** Necessary since gdb can't debug constructors. */
template <template<class> class chromatic_policy, 
	  typename T_rng_policy> 
void
mcrx::p04_grain_model<chromatic_policy, T_rng_policy>::
retarded (const std::string& data_directory, const Preferences& p,
	  const T_unit_map& units)
{
  // shortcuts
  std::vector<boost::shared_ptr<emitting_grain> >& grains = 
    grain_collection_grain_model<chromatic_policy, T_rng_policy>::grains_;
  std::vector<array_1>& dnda =
    grain_collection_grain_model<chromatic_policy, T_rng_policy>::dnda_;

  // copy units from those supplied into *scatterer* unit map.
  T_unit_map& u = 
    mcrx::scatterer<chromatic_policy, T_rng_policy>::units_;
  u["length"] = units.get("length");
  u["wavelength"] = units.get("wavelength");
  u["mass"] = units.get("mass");

  // The grains use their preferred units, which are also copied into
  // the grain model units.

  // Read the grains. The _v2 file is the actual cross sections from
  // the P04 file, which are slightly different from the ones
  // extracted from the Draine & Li files.
  grains.push_back
    (boost::shared_ptr<emitting_grain>
     (new thermal_equilibrium_grain
      (data_directory + "/p04dust_v2.fits", p)));

  T_unit_map& u2 = 
    mcrx::grain_model<chromatic_policy, T_rng_policy>::units_;
  u2["length"]=grains.front()->units().get("length");
  u2["wavelength"]=grains.front()->units().get("wavelength");
  // grains have no mass unit. But the size distributions do, so for
  // simplicity we select the scatterer mass unit as our mass unit.
  // NOTE that this arrangement means the size distributions and the
  // grains in general use different length units.
  u2["mass"] = u.get("mass");

  // we hardcode "size distribution" here, as it is just the number of
  // grains per mass.

  // this is the number of grains in 1Msun, ~7.6e46
  T_float npermass = constants::Msun/
    (4./3.*constants::pi*0.12e-6*0.12e-6*0.12e-6*3.6e3);
  // apply mass unit
  npermass *= units::convert(units.get("mass"),"Msun");

  // delta_size() of the grain will return 1. The "size distribution"
  // should be in scatterer units, so dnda here should be number of
  // grains per unit mass DIVIDED by the length unit conversion. This
  // will ensure that dn/da*da is n per unit mass.
  
  // evaluate the size distributions. Note that we need to convert
  // from grain length unit here to the scatterer units
  const T_float gsize_conv = units::convert(u2.get("length"),
					    u.get("length"));

  dnda.resize(grains.size());
  dnda[0].resize(grains[0]->asizes().size());
  dnda[0] = npermass/gsize_conv;


}

template <template<class> class chromatic_policy, 
	  typename T_rng_policy> 
mcrx::p04_grain_model<chromatic_policy, T_rng_policy>::
p04_grain_model (const std::string& data_directory, const Preferences& p,
		 const T_unit_map& units) :
  grain_collection_grain_model<chromatic_policy, T_rng_policy>()
{
  retarded(data_directory, p, units);
}


template <template<class> class chromatic_policy, typename T_rng_policy> 
class mcrx::p04_grain_model_creator_functor
{
public:
  typedef boost::shared_ptr<mcrx::grain_model<chromatic_policy, T_rng_policy> > return_type;
  
  p04_grain_model_creator_functor() {};
  return_type operator()(const Preferences& p) const {
    const std::string sigma_dir = 
      word_expand(p.getValue("grain_data_directory", std::string ()))[0];
    
    //THIS IS A HACK BECAUSE WE CAN ONLY HAVE ONE ARGUMENT
    T_unit_map units;
    units["mass"] = p.getValue("massunit", std::string());
    units["length"] = p.getValue("lengthunit", std::string());
    units["wavelength"] = p.getValue("wavelengthunit", std::string());

    return return_type
      ( new mcrx::p04_grain_model<chromatic_policy, T_rng_policy> 
	(sigma_dir, p, units));
  };
};


#endif
